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SUMMARY 


A method for solving the Navier-Stokes equations based on splitting the 
velocity vector into its rotational and irrotational parts has recently been 
applied successfully to internal flow computations. In this paper, the appli- 
cability of the method to external flows is examined by studying several model 
problems. The model problems are those of laminar and turbulent incompressible 
flow past a semi-infinite flat plate and laminar incompressible flow past a 
finite flat plate. For these problems, the procedure accurately reproduces the 
known solutions and is computationally very efficient even at high Reynolds 
numbers. Computational aspects of the method are discussed along with the 
possibility of using the procedure to retrofit a viscous capability into 
existing potential- flow codes. 


INTRODUCTION 

In recent years techniques for predicting the viscous flows over airfoils 
have been developed based on the numerical solution of the Navier-Stokes equa- 
tions (ref. 1). Although such solutions are useful, their speed of convergence 
is limited by the speed at which the basic numerical scheme can solve the Euler 
equations. For example, the method of Beam and Warming (ref. 2) solves the 
full Navier-Stokes equations in 20 percent more time than the same method takes 
for the Euler equations. An alternative to the direct solution of the Navier- 
Stokes equations has recently been proposed by Dodge (ref. 3) . According to 
Dodge, the velocity vector is split into its rotational and irrotational compo- 
nents. This split, coupled with an appropriate identification of the pressure, 
leads to a modified potential equation for the pressure field and a set of 
transport equations for the rotational component of the velocity. The poten- 
tial equation and the transport equations are then solved separately, with 
their coupling accounted for by iteration. The advantage of this method is 
that all the well-developed technology for solving the potential equation can 
be used for solving the pressure-field equation. Such a procedure can thus 
take advantage of the more rapid solution procedures available for the poten- 
tial equation compared with the Euler equations. 

The technique proposed by Dodge has primarily been applied to internal 
flows (ref. 3) . Consequently, an investigation was started to examine the 
suitability of the procedure for solving external flows with the goal of 
calculating viscous flow fields about airfoils. 

In the current investigation, several aspects of the procedure used by 
Dodge are investigated. The first question addressed in the current work is 
whether the procedure accurately recovers known solutions to the Navier-Stokes 
equations for laminar flows. Secondly, the suitability of the procedure for 
predicting turbulent flows is examined. Both of these questions are studied by 
calculating high Reynolds number flow over a semi-infinite flat plate. It is 
recognized that the semi-infinite plate problem can be handled well by conven- 



tional boundary-layer theory. The final phase of the stud/ is therefore 
involved with solution of a problem that cannot be treated by the usual 
boundary-layer methods. The problem of interest in this report is the 
detailed prediction of the laminar flow in the vicinity of the trailing edge 
of a finite flat plate. Computed results for all three problems studied are 
compared with known solutions. 

Use of trade names or names of manufacturers in this report does not con- 
stitute an official endorsement of such products or manufacturers, either 
expressed or implied, by the National Aeronautics and Space Administration. 


SYMBOLS 

A x constant in x-direction coordinate stretch 

Ay constant in y-direction coordinate stretch 

B coordinate scale factor (see eq. (13)) 

C constant in coordinate transformation 

Cf skin-friction coefficient 

F reciprocal of mesh interval in x-direction 

G reciprocal of mesh interval in y-direction 

I number of iterations to convergence 

M maximum mesh increment counter in x-direction 

N maximum mesh increment counter in y-direction 

P pressure coefficient with triple-deck scaling 

p pressure 

R Reynolds number based on free-stream conditions 

R x Reynolds number based on x 

r residual of continuity equation 

U x-component of total velocity vector with triple-deck scaling 

u x-component of total velocity vector 

u r x-component of rotational part of velocity vector 

V total velocity 


2 



V r rotational part of velocity 

v y-component of total velocity vector 

v r y-component of rotational part of velocity vector 

x longitudinal coordinate 

x-| longitudinal coordinate when mesh increment counter is 1 

x stretched x-coordinate 

x-j stretched x-coordinate when mesh increment counter is 1 

y normal coordinate 

y stretched y-coordinate 

a grid stretching parameter 

constants in free-stream potential function 
5 boundary-layer thickness 

e eddy viscosity 

T 1 normal coordinate in Gortler transformation 

X-j ratio of local wall shear to local Blasius theory wall shear 

K longitudinal coordinate in Gortler transformation 

p density 

$ irrotational velocity potential 

<j> perturbation irrotational velocity potential 

tj)^ unperturbed irrotational velocity potential 

X longitudinal coordinate with triple-deck scaling 

Subscripts: 

m mesh increment counter in x-direction 

n mesh increment counter in y-direction 

te trailing-edge location 

x,£ differentiation in x-direction 
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y,ri differentiation in y-direction 

00 free-stream conditions 
Superscripts: 

1 global iteration counter 

1 turbulent fluctuating quantity 

vector quantity 


PROBLEM FORMULATION 
Governing Equations 

In the present work the governing equations are the incompressible Navier- 
Stokes equations. These equations may be written as 

$ * V = 0 (1 ) 

(V • $)V = -^p - R -1 ^ x $ x V (2) 

where lengths have been nondimensionalized by a reference length L, velocity 

1 0 

by the free-stream velocity V , and pressure by — PV^ 2 . The velocity vector 

2 

V is split into rotational and ir rotational components according to the 
relation 

V = + V r (3) 

Further, the pressure is defined, as discussed in reference 3, in terms of the 
irrotational component as 


P 


= C - 


2 


(4) 


where C is an arbitrary constant. Substitution of equations (3) and (4) into 
equations (1) and (2) gives 

V 2 $ = • v r 
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(V r 


* V)V$ + _(fo> + v r ) • $]v r = -R~">$ X ^ x V r 


For turbulent flow in two dimensions, these equations can be time averaged 
and expanded to give 


^xx + $yy = “ (u r ) x " ( v r)' 


(5) 


u r$xx v r$xy (*^x ^ Uj-) (Uj-)j^ + ($y + Vr) (^iOy — R "^(^nJyy (V^$}jJ — (u'v 1 )' 


( 6 ) 


u r$xy v r$yy (^x ^r^ x ($y Vj.) (Vj-)y = R "^(Vr)yy + (V^$)y""| 


(7) 


The only component of the Reynolds stress tensor retained is the u'v* usually 
retained in boundary-layer theory. Streamwise diffusion is neglected. 

In the current work, with streamwise diffusion neglected, equations (6) 
and (7) are parabolic. However, equation (5), the pressure equation, still 
retains its elliptic character. Note that this system of equations has more 
generality than the usual boundary-layer equations; it admits a normal pressure 
gradient in the viscous layer. It is also more general than the usual "parabo- 
lized" Navier-Stokes equations; it does not require a split of the pressure 
into a known marching direction pressure gradient and a viscous perturbation 
pressure field. (See, e.g., ref. 4.) 


Boundary Conditions 

As indicated in the Introduction, the problem of interest in the current 
investigation is laminar and turbulent flow over semi- inf inite and finite flat 
plates. The computational domain considered and the associated boundary condi- 
tions on the primitive variables are shown in figure 1 . The leading edge of 
the plate was excluded from the computational domain. Therefore, techniques 
for dealing with the leading-edge singularity did not have to be considered 
(ref. 5) . This subject is left to a future investigation. The inflow boundary 
of the computational domain is taken sufficiently far downstream of the plate 
leading edge that effects of the plate leading-edge singularity have vanished 
and the flow is well described by the Blasius boundary-layer profile with an 
undisturbed external flow (which is not precisely accurate, but is accurate to 
R x _1 /2 f or large enough distances from the leading edge) . 
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The boundary conditions on $ and V r are obtained from the boundary 
conditions on the primitive variables as follows; 

Wall boundary 



Wake centerline boundary 


(u r ) y - v r 
<5 y = 0 



(y = o, x s x te ) 


(y = 0, x > x te ) 


Since the flow is irrotational outside the boundary layer, the rotational com- 
ponent of velocity must vanish outside the boundary layer. 

y -> <*> boundary 


u r = v. 


$ 


= 4), 



(y **■ °°r H X ^ 


where ^ represents some as yet unspecified potential function. 
x ■> oo boundary 


$ = cjj^ (x *► °°, 0 ^ y ^ °°) 

No boundary conditions are required on u r and v r as x -*■ 00 , since the 
transport equations governing their evolution are parabolic. 

The boundary conditions along the line x = 1 pose a particular problem 
for the present method, since the splitting of the known total velocity vector 
into rotational and irrotational components is not unique. The initial choice 
made in the present investigation was the Blasius velocity profile and zero 
longitudinal pressure gradient. 

Inflow boundary condition 


(x = 1, 0 = y 
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u r = u ” (<Ux 


(8a) 



V r = V - (Cjjjy 

II 

— 1 

o 

HA 

HA 

(8b) 

o 

II 

(x = 1, 0 ^ y ^ 

(8c) 


However, computational results forced modification to these conditions which 
are discussed in some detail in a subsequent section. 


Turbulence Model 

The Reynolds stress term in equation (6) is modeled with the usual eddy 
viscosity assumption; that is, 


3u 

-u'v* = £ — 
3y 


Substitution of the x-component of equation (3) into the preceding equation 
gives 

-u'v' = e(u r ) y + (9) 

Substitution of equation (9) into equation (6) leaves the specification of e 
to close the system of governing equations. In the present work, a two- layer 
eddy viscosity formulation is used, as discussed by Cebeci and Smith (ref. 6) . 
In the inner region the formulation is based on Prandtl's mixing- length model 
along with the Van Driest damping factor. In the outer region, Clauser's 
velocity defect model is used. No normal intermittency is used. 


SOLUTION OF GOVERNING EQUATIONS 
General Solution Strategy 

In the present work, the governing equations, equations (5) to (7), are 
solved with a finite-difference procedure. The general strategy of the 
solution, as proposed in reference 3, is to solve equation (5) separately from 
equations (6) and (7) , and to account for their coupling iteratively. This 
procedure allows the use of the technology available for solving the potential 
equation, thus treating the pres sure- field solution as an elliptic problem, 
while still solving the viscous portion of the flow field with a marching 
calculation. 

In practice, the solution procedure is started by first determining the 
inviscid pressure field by specifying an initial potential function which sat- 
isfies equation (5) with the right-hand side set to zero. Once the $ field is 
established, the u r and v r fields are established by solving equations (6) 
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and (7) with a marching solution initiated by the assumed profile at x = 1 . 
All the required 0 derivatives in equations (6) and (7) are known functions 
and are calculated from^the solution of equation (5) . Once the u r and v r 
fields are found, V • V r can be calculated and treated as a known source 
term for the solution of equation (5) . This process is continued until the 
continuity equation is satisfied throughout the field to some required degree 
of tolerance; that is, 

| r 1 1 = | (V2$)I + (V • V r ) 1 [ < Tolerance (1 < x < 00 , 0 ^ y < °°) (10) 

where the superscript I refers to iteration number. 


Transformation of Momentum Equations 

Examination of equations (6) and (7) reveals that they need only be solved 
in regions where the velocity has a nonzero rotational part; that is, in the 
viscous layer. In most situations the viscous layer has the characteristic 
that its thickness increases in the streamwise direction. For this reason, the 
boundary-layer equations are usually transformed into a coordinate system in 
which the height of the viscous layer is approximately constant. For compara- 
tive purposes, a similarity transformation is introduced in the current inves- 
tigation by introducing Gortler scaling on the normal coordinate. 

The x-coordinate is now alined with the free-stream direction and the fol- 
lowing new independent variables are introduced: 


K = x 


n 



where 


B(5) = 1 


( 11 ) 


( 12 ) 


(13a) 


if the normal coordinate is not transformed, and 


B(5) = (2?)V2 


(13b) 


if the normal coordinate is transformed by Gortler scaling. In these new 
variables, equations (6) and (7) become, after substitution of equation (9), 
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s2 ( u r $ xx + v r $ xy) + B 2 (u r + $x)( u r)£ + [ B < v r + $ y) “ Bn(u r + $ x )j (u r ) 


= R' 


'^[o + e) (« r )n] n + B 2 (V 2 $) X - B 2 (£$ xy )y| 


( 14 ) 


B 2 (u r 4 , xy + v r $yy) + B 2 (u r + <J> X ) (V r )£ + J^B(v r + 4>y) - Br|(u r + 4> x )j (v r ) n 


= R-l[(v r ) nn + B 2 (V2$) y J (15) 


where B = (B - 1 ) / (2£ - 1 ) . 


Numerical Method 

Solutions to the system of equations (5), (14), and (15), with their asso- 
ciated boundary conditions, are obtained with finite-difference methods. Since 
equation (5) and equations (14) and (15) are of different types, the numerical 
techniques used in their solution differ. Equation (5) is solved by the 
successive line overrelaxation method (SLOR) discussed in reference 7; equa- 
tions (14) and (15) are solved by using an implicit marching technique. 

Consider first the finite-difference mesh used in the solution. In all 
cases, the equations are solved on nonuniform meshes. Stretched Cartesian 
coordinates are used in the solution of the potential equation (eq. (5)) for 
all cases, and are used in the solution of the viscous equations (eqs. (14) 
and (15)) when B = 1. The general mesh notation is indicated in figure 2. 

The nonuniform Cartesian mesh is obtained by coordinate transformations of 
the form 


A x x 



(16) 


Ay y 

y = (17) 

1 - y 2 


where Xq = 0 for the semi-infinite flat-plate case and x Q = xt e for the 
finite flat-plate case. The actual computations are carried out in the x,y 
plane with a uniform mesh. The step sizes in the x- and y-directions are 
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Ay = 



where C = 1 for the semi-infinite plate, C 


plate, and x-j = x^e - 1. 


+ 


X 1 


for the finite 



+ x-| 2 


When B = (25) 1 / 2 , different meshes are used in the solution of equa- 
tion (5) and equations (14) and (15). For this case, the stretched Cartesian 
grid is used for equation (5), whereas the n-coordinate in equations (14) 
and (15) is stretched using the Roberts transformation (ref. 8). This trans- 
formation can be written as 



where a is a stretching parameter. Equation (16) is still used to stretch 
the x-coordinate. 

New consider the solution to equation (5) . In the present work, equa- 
tion (5) is solved for a disturbance potential by writing 


V 2 $ = ^ 2 (0 OO + <t>) = -V • V r 
and, since 6 is chosen such that 

CO 

V 2 ^ = 0 

then 

V 2 <)> = • v r 


(18) 
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Now Fjjj and G n are defined as 


■m 


(Ax) - 1 


dx 

dx 


G n = (Ay )" 1 


dy 

dy 


With this notation, equation (18) is discretized according to 


^xx [(^m+ljn “ ^n^m+O/^) “ ( < J ) m,n " ^m-1 ,n) F m- (1/2) F m ( 19a ) 


^yy = £^m,n+l " ^m,n) G n+(l/2) “ (§m,n “ ^in,n-l > G n- (1/2) ] G n 09b) 


The boundary conditions are discretized according to 


'I’m,! 
h ,n 
^m,N 
^M,n 


3 

^2,n 

0 


= 0 


(2 ^ m ^ M - 1 ) 
(1 ^ n ^ N) 
(1 ^ m ^ M - 1) 
(2 S n ^ N) 


The discretized representation of equation (18) is solved by the SLOR procedure 
with y as the implicit direction and with the procedure sweeping from n = 2 
to n = N - 1 on each iteration cycle. The most recent values of 4> are used 
whenever possible. The Thomas algorithm is used to solve the set of implicit 
tridiagonal equations along each line where m is constant. If r m n is 
defined as ' 


I (^^mfn + $ * v r)m,n 
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then convergence is achieved when 


r m,n 


$ — (Ax ^ + Ay^) 
2 


(2 ^ m ^ M - 1 , 2 = n = N - 1 ) 


Now consider the solution of the rotational-velocity-component momentum 
transport equations, equations (14) and (15). These equations are parabolic 
in the £ -direction and hence are solved by downstream marching from an ini- 
tial profile located at £ = 1 . If g denotes either u r or v r , then 
equations (14) and (15) are discretized according to 


1 

(9x)m,n = ~^9m,n " ^9m-l ,n + 9m-2,n)*'m (20a) 

1 

(9y)m,n = T^9m f n+1 ” 9m,n-l) G n (20b) 

(9yy)m,n = £(9m,n+l “ 9m, n) G n+ (1/2) “ (9m, n ” 9m,n-l) G n-(l/2) J G n (20c) 


If these equations were substituted directly into equations (14) and (15), a 
set of N - 2 simultaneous nonlinear algrebraic equations would result along 
the mth column. Thus, before discretization, equations (14) and (15) are 
quasi-linearized. When f is either a u r or v r derivative, then 


(^9)m,n = (^9)m,n + (^9^m,n “ (^9)m,n 


where the tildes mean guessed values. Substitution of these relations into 
equations (14) and (15) yields a set of 2 (N - 2) simultaneous linear alge- 
braic relations. The equation set can be completed by using the boundary 
conditions 


( u r)m,l ~ " (^*x)m,l 


( v r)m,l - - (^y)m,l 
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< u r)m,N “ 0 


( v r>m,N = 0 

This system of equations is then solved by the Gaussian elimination for a 
2x2 block tridiagonal matrix. 

Quasi-linearization of the viscous term requires special attention in 
turbulent cases. The nonlinear viscous term in equation (14) can be written 
as 



where a (5) is a function of boundary-layer thickness 6 and b is a con- 
stant (ref. 9) . The quasi-linear form is then 



3u r 

As can be seen, this form does not properly quasi-linearize the term a (6) , 

3ri 

since <5 is a function of the local solution. With this form, quadratic con- 
vergence is obtained on the column iterations for the laminar cases, whereas 
for the turbulent cases a 20-percent decrease in computer time over Picard 
iteration is achieved. 

In the solution of equation (5) , V • v r is treated as a known quantity 
at each mesh point. In the solution of equations (14) and (15), however, 
and all its required derivatives are treated as known. The quantity V • V r 
is numerically evaluated for each rotational-velocity-coraponent calculation by 
using equations (20a) and (20b) . The quantities $ xx and $yy are calculated 
after each solution for $ by using equations (19) along with known values of 
an ^ ^°cyy* In addition, the following relationships are used to obtain 

other required derivatives of 4>. 


(^x)m,n ~ “ < f > m-l,n) F ra + 
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( $ y)ra,n “ 2 ('t’m, n+1 “ 4>m,n-l) G n + (^ooy) rafn 
(^xy)m,n = “(^m+l ,n+1 ” ^m— 1 ,n+1 + ^m-l ,n-1 ~ ^m+l ,n-l ) F m G n + ^cojjyj mfn 

(V 2 $) m#n = (^xx^jii + ^yy^niin 
[(V 2 *)J = "rt ^ 2( t , )m+l ,n “ (^^m-l ,nJ F m 

£(V 2< $) y J = tE (^ 2 4>) m,n+1 “ (^ 2< t>) m,n-lJ G n 

u -Tn, n 2 


Since these quantities are required only on interior mesh points, special 
formulas on the boundaries are not required. 

It should be noted that in the actual programming of the scheme, provision 
was made for resolving the rotational component of velocity on a finer mesh than 
the irrotational component. When B = 1 , this was done by breaking up the x- 
and/or y-increments in the computational plane into the desired integer number 
of subincrements. When this provision was made, the derivatives of $ required 
in the solution of equations (14) and (15) were linearly interpolated onto the 
finer grid. 

When Gortler variables are used for solution of equations (14) and (15), 
linear interpolation on both 0 and V ♦ v r is used. 


RESULTS AND DISCUSSION 

As discussed previously, three test cases were chosen for computation with 
the split Navier-Stokes formulation. These cases are the laminar and turbulent 
semi-infinite flat plate and the finite flat plate. The computer results for 
these test cases are discussed separately in the following sections. 


Laminar Semi-Infinite Plat Plate 

It is well known that a singularity exists in the pressure and vorticity 
at the leading edge of a flat plate (ref. 5) . Rather than deal with this 
problem in the current investigation, all solutions were initiated a suffi- 
cient distance downstream of the leading edge such that the effects of this 
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singularity were negligible. According to reference 5, at a Reynolds number 
based on distance from the leading edge R x of 10 4 the boundary layer is well 
described by the Blasius velocity profile. For the current computations, the 
free-stream conditions were such that R x = 10 4 for x = 1. The initial com- 
putations were started at x = 1 with the Blasius velocity profile. This pro- 
file was generated in the program itself by numerical solutions of the appro- 
priate nonlinear ordinary differential equations (ref. 10). 

In order to initiate the calculation, a choice had to be made for <b , the 
basic potential function. The most obvious choice to make was 


<b = x + Constant 

' CO 


( 21 ) 


since this represents the undisturbed inviscid flow about a flat plate. This 
choice led to immediate difficulties in equation (8b) since (^y is zero 
everywhere. Use of the Blasius-predicted v in equation (8b) would then lead 
to nonzero values of v as y ■+ 00 , which is physically incorrect. For this 
reason, v was arbitrarily set to zero along the initial profile when this 
form was chosen for d> . 

‘ OO 

A second possible choice for ^ is that for flow about the displacement 
body presented by the plate boundary layer to the inviscid flow. It is well 
known that the flat-plate displacement body is a parabola which has infinite 
thickness as x °°. In fact, the velocity potential is infinite as y °° 
for such bodies, whereas the boundary conditions and equation (21) indicate a 
finite value for the potential function as y °°. One way to alleviate this 
difficulty is to choose the potential function for flow about the Blasius 
displacement-body parabola for (j> ; that is, 

BR-l/2 

4> = x + 

1/5 T 2 2 V2l V2 

|_x + (x 2 + y 2 ) ' J 


where 3 = B b = 1.72078765 for the Blasius displacement body (ref. 11). 
Note that when 3=0 the uniform-stream potential function is recovered. 
The preceding potential function predicts nonzero values of (^Jy along 
the initial data line; that is, 


^co)y 


BR-V 2 [x + 
2\/2 


— jl / 2 

(x 2 + y 2 )^/ 2 
(x 2 + y 2 ) 


15 



Thus, the normal- velocity profile on the inflow boundary is assumed to have the 
following form: 



Normal— velocity components 


v r = v Blasius “ y v Blasius - ^oo^y 

v r = v Blasius > (^oo)y 

Also, since (^y ¥■ 0 along y = 0, the surface boundary condition becomes 

v r = -<<U y (Y = 0, 1 S x < «) 


Calculations were made for the semi- inf ini te plate using these two 
free-stream potential functions on a mesh with 41 points in the x-direction 
and 58 points in the y-direction, and with A x = 20 and Ay = 0.22. For 
the case where 3 = 3g, the calculation was stopped when the convergence 
criterion was satisfied, a total of 44 global iterations (i.e., iterations 
between eq. (5) and eqs. (14) and (15)). Convergence was extremely slow for 
the 3=0 case, however ,^.and the calculations were arbitrarily stopped after 
200 iterations with |V • v| max = 6.48 x 10“ 4 and the wall shear converged 
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to four decimal places of accuracy. The distribution of X-| , the ratio of pre- 
dicted wall shear to the local Blasius value, along the plate is presented in 
figure 3 for these cases. As can be seen, the predicted wall shear immediately 
departs from the Blasius value in the vicinity of x = 1 in both solutions. 
This is attributed to the arbitrary split of the initial velocity distribution 
into its rotational and irrotational parts. For the case with 3 = B b , the 
Blasius wall shear is quickly recovered further down the plate, whereas the 
case with 3=0 never recovers the Blasius solution. This behavior is 
attributed to the fact that the solution obtained with 3 = 3 B accounts for 
the interaction, however weak, between the boundary layer and inviscid flow. 

The choice of 3 = 0 coupled with the boundary condition as y -* 00 ignores 
this interaction and forces the potential function to vanish as y °°. This 
leads to serious errors, as can be seen in figure 4, in the prediction of the 
normal-velocity component all the way down to the wall. The conclusion drawn 
from these calculations was that for interaction calculations involving infi- 
nite bodies, the proper choice of boundary conditions at y °° is crucial. 

No further calculations were made with 3=0 for the semi-infinite plate. 

After establishing the ability of the split equations to accurately 
reproduce the Blasius solution igith coincident meshes for both the pressure 
(i.e., 4>) and viscous (i.e., V r ) solutions, calculations were performed with 

increasingly coarsened (in the y-direction) pressure meshes while holding the 
viscous mesh constant. This was accomplished by holding A y the same on both 
grids and successively halving the number of points in the y-direction for the 
pressure calculation. The resulting pressure grid points all coincided with 
viscous grid points with an integral number of viscous grid points between. 
Three successive halvings of the number of pressure-mesh points were done over 
the 41 x 58 mesh. The results of these calculations are presented in table 1 
in terms of the percent error in friction drag relative to the Blasius value 
over the interval 1 = x = 90.12, along with the number of iterations to con- 
vergence and computer time required. All runs were converged to the truncation 
error of the pressure mesh. Little accuracy is lost by the coarsening of the 
pressure mesh, and there is a substantial savings in computer time. The con- 
vergence histories of these calculations are presented in figure 5. Coarsening 
of the pressure mesh is shown to accelerate the overall convergence rate. 

To further assess the potential of the velocity-split formulation, addi- 
tional calculations were carried out with entirely independent meshes for the 
pressure and viscous calculations by choosing B = (2?) '/ 2 . This allowed the 
number of viscous mesh points in the boundary layer to be held constant at 
about 40 points over the entire length of the plate. In this case, the viscous 
grid was again held fixed and A y was varied. Results of these calculations 
are presented in table 2. In this table, the ratio of the pressure-mesh incre- 
ment to the viscous mesh increment at the wall is used as a measure of the grid 
coarsening. It should be noted that in the solutions for mesh ratios of 32.90 
and 65.79 (A v = 3.52 and 7.04) the global iteration had to be under relaxed in 
order to achieve convergence. The reason for the deterioration of the solution 
accuracy on the coarser pressure grids is obvious from figure 6. In this fig- 
ure, the distribution of V • v r across the boundary layer is plotted 
along with the location of both the viscous and pressure-mesh points. As can 
be seen, when the pressure grid i§ coarsened to the point that it can no longer 
resolve the distribution of V • v r , serious deterioration in accuracy 
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occurs. Also presented in table 2 are the error in the local wall shear at 
x = 33.33 and the number of pressure-mesh points in the boundary layer. These 
data can be interpreted as a first-order estimate of the number of pressure 
mesh points that are required in the boundary layer to give adequate resolu- 
tion. From table 2 it can be seen that between 5 and 10 pressure-mesh points 
are required within the boundary layer to accurately predict the local wall 
shear. 


Turbulent Semi-Infinite Flat Plate 

The semi- infinite plate was also run for the turbulent case. In this 
case, the value of R x was 1.0 x 10 6 at x = 1.0. The solution was started 
from the Blasius solution at x = 1.0, and instantaneous transition was arbi- 
trarily introduced at x = 2.44 (m = 4). For this case the value of 0 was 
the laminar value since the exact potential- flow solution about the turbulent 
displacement body was not known. 

The distribution of skin friction along the plate is compared in figure 7 
with the Prandtl-Schlichting formula for a smooth flat plate (ref. 10). For 
the cases shown the grid was defined by the following parameters: 


Viscous 


Inviscid 


M = 41 


M = 41 


N = 54 


N = 61 


a = 2048 
r| e = 0.08 


Ax = 20 


Ay = 0.2, 0.3, 0.4 
A x = 20 


The solutions for the three values of Ay cannot be differentiated on 
the scale of figure 7. For values of Ay > 0.4, however, the solutions 
diverged while values of Ay < 0.2 left too few points outside the boundary 
layer on the inviscid mesh, and hence were not run. The three cases of 
Ay = 0.2, 0.3, and 0.4 took 17, 15, and 54 global iterations requiring 35, 

30, and 102 seconds, respectively, on the Control Data CYBER 175 computer to 
achieve convergence. 

The turbulent case thus appears somewhat more sensitive to the resolution 
of the boundary layer by the pressure calculation than the laminar case. For 
the case of Ay = 0.4, the first pressure grid point off the wall corresponds 
to y/6 =0.37 at x = 2.44, the first turbulent profile. For the turbulent 
case, resolution coarser than this leads to divergent solutions rather than 
inaccurate solutions as in the laminar case. On the other hand, these results 
show that the pressure grid does not have to resolve the details of the turbu- 
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lent boundary layer near the wall in order to produce accurate results, pro- 
vided the viscous mesh is fine enough. 

In view of the results for the laminar case, the apparent accuracy of the 
turbulent calculations is surprising since the basic potential function does 
not properly account for the decay of the normal component of velocity as 
y °°* It should first be noted, however, that the predicted solution in 
figure 7 does begin to depart from the data curve fit at the higher values of 
R x and examination of the normal-velocity profiles in this vicinity does show 
some small negative values near the wall. Secondly, the values predicted by 
any finite-difference program of this type depend upon the transition history 
and turbulence model employed. Thus, for this case an exact solution is not 
available for comparison, as in the laminar case. Hence, a direct quantitative 
error estimate is not possible. In addition, it can be seen from figure 7 that 
the predicted values of skin friction are in error in the laminar region. This 
error is attributed to the upstream influence of the transition region through 
the pressure calculation. 

It should finally be noted that for turbulent calculations the Newton 
iteration in the viscous marching-column solution did not converge quadrati- 
cally. This is in agreement with the results of reference 9. 


Finite Flat Plate 

The third problem chosen for solution in the current investigation is the 
laminar finite flat plate for R — 5 x 10^. In this problem, a singularity 
exists at the trailing edge of the plate due to the discontinuous boundary con- 
dition. Calculations were made with the present method using coincident pres- 
sure and viscous meshes with 3=0. Results of these calculations are shown 
in figures 8 to 10. In these figures, the following triple— deck scaled quanti- 
ties are presented: 


X 5 /4( X _ x te ) 

X = 

e3jc te 


u 

U = 

eXl/4 

P - Pa, 
P = 

- e2 X l/2 
2 


e = r~1/8 


(X = 0.33206) 
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These results are compared with the triple-deck theory solutions of Melnik 
and Chow (ref. 12) and Jobe and Burggraf (ref. 13). Figure 8 presents the 
wall/centerline pressure distribution in the vicinity of the trailing edge, 
figure 9 the wall shear, and figure 10 the wake centerline velocity distribu- 
tion. Agreement between the present calculations and triple-deck theory is 
excellent. The better agreement between the present results and those of 
reference 13 is attributed to the fact that the results of reference 12 were 
calculated on a much finer grid than that used in either the current calcula- 
tions or those of reference 13. Additional calculations were made with the 
pressure solution obtained on every second and fourth viscous mesh point with 
no plottable difference in the results. The calculations required 75 seconds 
of CYBER 175 time for the coincident mesh case, whereas the two coarser grid 
calculations required 50 and 39 seconds, respectively. 


CONCLUDING REMARKS 

It is concluded from the current investigation that the split- velocity 
Navier-Stokes formulation can be used to obtain accurate solutions to the 
Navier-Stokes equations for finite and semi-infinite flat plates. Results of 
the pressure-mesh coarsening studies carried out in this investigation also 
indicate that the method can be used to retrofit a viscous capability into 
existing potential-flow codes. In order to retrofit such a capability, how- 
ever, care must be taken to adequately resolve the viscous layer with the pres- 
sure grid. Since existing potential-flow computer programs have computational 
grids selected to resolve purely inviscid phenomena, only in special cases can 
the retrofit be made without altering the grid. For the commonly used numeri- 
cal schemes in modern potential- flow programs, the introduction of a known 
function on the right-hand side of the potential equation poses no difficulty. 

For the problems studied in the current investigation, the computational 
efficiency of the method appears unaffected by Reynolds number. The reason 
is that in any given calculation the bulk of the computational effort (about 
80 percent) is in the viscous marching calculation. In this portion of the 
calculation, the effect of Reynolds number is accounted for by scaling the 
viscous grid. Hence, the total computation time for a given viscous marching 
sweep is unaffected by Reynolds number. Finer resolution of the viscous region 
by the pressure grid at higher Reynolds number does slow this portion of the 
calculation down but has little effect on the overall computation time. For 
this reason, the method appears promising for calculating steady, high Reynolds 
number, external flows. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
March 18, 1980 
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TABLE 1 EFFECT OF PRESSURE-MESH COARSENING ON 


OVERALL SOLUTION ACCURACY FOR RELATED MESHES 


Number of 
y- points 
(pressure) 

Percent error 
in friction 
drag 

I 

CYBER 175 
time, 
sec 

58 

O 

in 

• 

o 

45 

47 

30 

o 

• 

T 

25 

20 

16 

3.70 

12 

8 

9 

3.73 

9 

6 


TABLE 2.- EFFECT OF PRESSURE-MESH COARSENING ON 
SOLUTION ACCURACY FOR INDEPENDENT MESHES 


Wall mesh 
ratio 

Percent error in 
friction drag 

x = 33.33 

Number of y-points in 
boundary layer (pressure) 

X 1 

2.06 

-0.49 

46 

0.998 

4.11 

-.70 

32 

.996 

8.23 

-.16 

19 

.998 

16.45 

2.46 

10 

1.029 

32.90 

13.4 

5 

1.142 

65.79 

23.0 

3 

1.327 
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(a) Semi-infinite flat plate. 

Figure 1 Coordinate system and boundary conditions for test problem. 













Figure 5.- Effect of pressure-mesh coarsening on 
convergence of global iteration. 
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